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General relativistic superfluid neutron stars have a significantly more intricate dynamics than 
their ordinary fluid counterparts. Superfluidity allows different superfluid (and superconducting) 
species of particles to have independent fluid flows, a consequence of which is that the fluid equations 
of motion contain as many fluid element velocities as superfluid species. Whenever the particles 
of one superfluid interact with those of another, the momentum of each superfluid will be a linear 
combination of both superfluid velocities. This leads to the so-called entrainment effect whereby 
the motion of one superfluid will induce a momentum in the other superfluid. We have constructed 
a fully relativistic model for entrainment between superfluid neutrons and superconducting protons 
using a relativistic a ~ u> mean field model for the nucleons and their interactions. In this context 
there are two notions of "relativistic" ; relativistic motion of the individual nucleons with respect to 
a local region of the star (i.e. a fiuid element containing, say, an Avogadro's number of particles), 
and the motion of fluid elements with respect to the rest of the star. While it is the case that the 
fluid elements will typically maintain average speeds at a fraction of that of light, the supranuclear 
densities in the core of a neutron star can make the nucleons themselves have quite high average 
speeds within each fluid element. The formalism is applied to the problem of slowly-rotating 
superfluid neutron star configurations, a distinguishing characteristic being that the neutrons can 
rotate at a rate different from that of the protons. 



I. INTRODUCTION 



A new generation of gravitational wave detectors (LIGO, VIRGO, etc) are now working to detect gravitational 
waves from compact objects, such as black holes and neutron stars. With this detection we expect to have a unique 
probe of the physics that dictates their behavior. This is ushering in a new era where strong-field relativistic effects 
will play an increasingly important role. Only through their inclusion can we hope to accurately decipher what 
gravitational wave data will have to tell us. With that in mind, we present here a fully relativistic model of the 
so-called entrainment effect (to be described in some detail below) that is a necessary feature of the dynamics of 
superfluid neutron stars. 

For the densities appropriate to neutron stars there are attractive components of the strong force that should 
lead, via BCS-like mechanisms, to nucleon superfluidity and superconductivity. Indeed, calculations of supra-nulcear 
gap energies consistently lead to the conclusion that superfluid neutrons should form in the inner crust of a mature 
neutron star, with superfluid neutrons and superconducting protons in the core. Even more exotic possibilities have 
been suggested, such as pion condensates, superfluid hyperons, and superconducting quark matter. Perhaps most 
important is the well-established glitch phenomenon in pulsars the best description of which is based on superfluidity 
and quantized vortices. Superfluidity should affect gravitational waves from neutron stars by modifying the rotational 
equilibria and the modes of oscillations that these objects support [0, |[ |j. 

The success of superfluidity in describing the glitch phenomena is due in part to the fact that the superfluid 
neutrons of the inner crust represent a component that can move freely (for certain timescales) from the rest of the 
star. Explaining the glitch phenomena then becomes a question of how to transfer angular momentum between 
the various "rotationally decoupled" components. For the modes of oscillation, it is by now well established that a 
similar "decoupling," this time between the superfluid neutrons of the inner crust and core and a conglomerate of the 
remaining charged constituents (e.g. crust nuclei, core superconducting protons, and crust and core electrons), leads 
to a mode spectrum for superfluid neutron stars that is quite different from that of their ordinary fluid counterparts 
(see [Q, and references therein, for a complete review). 

Several recent studies H, ||, |j have established that the entrainment effect is an important element in modelling 
the rotational equilibria and modes of oscillation of superfluid neutron stars. Sauls |^ describes the entrainment 
effect as a result of the quasiparticle nature of the excitation spectrum of the superfluid and superconducting nucleons. 
That is, the bare neutrons (or protons) are accompanied by a polarization cloud containing both neutrons and protons. 
Since both types of nucleon contribute to the cloud the momentum of the neutrons is modified so that it is a linear 
combination of both the neutron and proton particle number density currents, and similarly for the proton momentum. 
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Thus when one species of nucleon acquires momentum, both types of nucleons will begin to flow. 

In the core of a neutron star, the Fermi energies of nucleons (as well as some of the leptons) can become comparable 
to their mass-energies, because the Fermi energies are a function of the local particle number densities, and these 
can be quite high. This implies that any Newtonian model for entrainment must become less reliable as one probes 
deeper into the core of a neutron star, and thus a relativistic formulation is required. In fact, we will see that the 
Newtonian parameterized model of Prix et al. does deviate most from the relativistic model in the core. There 
are two purposes for which a relativistic formulation is necessary. At the microscopic level, the nucleons will (locally) 
have average speeds that are comparable to the speed of light. As well, at a mesoscopic level, the fluid elements, 
which contain a large number of nucleons, could have average speeds that are also comparable to the speed of light. 
The formalism that we develop here will be relativistic in both respects. One should note, though, that in realistic 
astrophysical scenarios (e.g. when an isolated neutron star undergoes linearized oscillations, or a pulsar exhibits a 
glitch) the fluid element average speeds are typically only a few percent of that of light. 

To date, studies of superfluid dynamics in neutron stars have relied on models of entrainment that are obtained in 
the Newtonian regime. For instance, a few of the most recent studies |l^ have employed a parameterized model 
for entrainment that is inspired by the Newtonian, Fermi-liquid calculations of Borumand et al. An alternative 

formulation — motivated by mathematical simplicity that allows for analytic solutions for slowly rotating Newtonian 
superfluid neutron stars — for parameterizing entrainment has been recently put forward. Here we take a different 
approach, and this is to use a ct — w relativistic mean field model, of the type that is described in detail by Glendenning 
[g3[. Although a relativistic Fermi- liquid formalism exists |Q, we prefer to use the mean field model because it is 
sufficiently simple that semi- analytical formulas result, and a clear connection between the coupling parameters at 
the microscopic level can be made to the macroscopic properties (such as mass and radius) of the star. An immediate 
consequence is the ability to compare the relativistic entrainment model with the two parameterized models. We will 
see that the model used by Prix et al. although limited, is a better fit than the other formulation. 

The next section begins with a review of the a — uj model. That is followed by an application of the mean field 
approximation to obtain an equation of state that includes entrainment. In Sec. 3, we briefly review the general 
relativistic superfluid formalism and how it is used to describe slowly rotating configurations. We then use the mean 
field results to produce explicit models. After some concluding remarks, an appendix is given that contains some of 
the technical details and results. Throughout we will use "MTW" |jl^ conventions, a consequence of which is that 
several equations will have minus sign differences with, for instance, those of p3|. 



II. RELATIVISTIC MEAN FIELD THEORY OF COUPLED FLUIDS 



To create a seamless conceptual basis for general relativistic calculations of dynamic processes in neutron stars, we 
need a covariant formalism that describes the strongly interacting coupled neutron and proton fluids. It should be 
sufficiently simple that it provides physical insight, yet accurate enough that it can serve as the basis for realistic 
numerical calculations. For static stars, this role is played by the a — lo effective mean-field theory [Q. Our task 
in this paper is to generalize this theory to dynamic stars. In particular, we are interested in situations where there 
is relative motion of the two fluids, since the entrainment of one by the other turns out to play a large role in the 
dynamics. 

The Lagrangian density for the baryons and the mesons that the baryons exchange is as in the static case. It is 

L = Lb + La + + Lint , (1) 

with 



Lb - ^Jiij^d'' - m)4, (2) 

as the baryon Lagrangian. Here ip is an 8-component spinor with the proton components as the top 4 and the neutron 
components as the bottom 4. The 7^ are the corresponding 8x8 block diagonal Dirac matrices. The Lagrangian 
for the CT mesons is 



1„ 1 



The Lagrangian for the lo mesons is 



La = --d.ad^a - -mia' . (3) 



Lu^ = -iw^^w'^'' - im^w^w^ (4) 
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where co^v = d/^uii, — d^ui^. The interaction Lagrangian density is 

Lint = ga(r'ipil> - guj'^n'4>i^'4} ■ (5) 

The Euler-Lagrange equations are 

(-□ + ml)a = g^i,^) , (6) 

(-□ + ml) + d^^d'^LOv = -gu^i'lu'tp , (7) 

ihud'^ - m)tlj = goj^nUJ^i) - QaC^"^ ■ (8) 
Finally, the stress-energy tensor takes the form 

containing contributions from the baryons (6), the mesons {a,uj), and the interaction. Individually, these are 

T^"" = -iV;(7^a'' - r]^"'j"dc)^ - 1717}^"^^) , (10) 

T^-' = d^'ad'^a - iry^-^m^^ - ^r]^"'d"'adaa , (11) 

= id^uj'^ - d'^uj^) 5^a;„ - ^^''m^a;"^;^ - ^v^^mlu^^uj^p , (12) 

TZ = v^^gMi^ - n'^^'g.ujo.i^ri, . (13) 

We now solve these equations in the mean field approximation, eventually in a frame in which the neutrons have 
zero spatial momentum while the protons have on average a wavevcctor K^^ = (Ko,0,0, Kz). In this approximation 
we ignore all gradients of the averaged sigma and omega fields, and the neutrons and protons are taken to be in 
plane-wave states. The problem simplifies considerably and we find for the a and fields and the stress-energy 
tensor Tjf that 

m* = m-cl {'ipi)) , (14) 

{gu;Oj^) = -cl{^j^^) , (15) 

= -I (cj2 {g^uj'^) {g^cj^) + [m - m.f) Sj^ - i (*7''5.*) , (16) 

where, for later convenience, we have introduced the notation = {ga/mfrY a^^d = (ffw/^^w)^ and the Dirac 
effective mass m*, i.e. 

{gcru) =m-m^ . (17) 
Restricting to the zero-momentum frame of the neutrons leads to a set of algebraic equations for the field: 

(ffa;Wo) = cK^-iolli) , (18) 

{gu^^z) = cl{ip^,i>) . (19) 

The final equation is not needed in the case where both neutrons and protons have zero average momentum, since (wz) 
then vanishes by isotropy. In this case, the neutrons and protons have a common rest frame and ('tp^^ip) = tp^tp = n+p 
where n and p are the baryon number densities of the neutrons and protons, respectively. The addition of the spatial 

velocity component complicates the solution of the problem considerably, in part because there is no longer a common 
rest frame for all the baryons. Each expectation value on the RHS of these equations involves an integration over 

1 /3 

the Fermi spheres of the particles, whose radii can be shown (c.f. the next section) to be /c„ = (37r^n°) and 

1 /3 

kp = (Stt^p") , where nP (p°) is the zero-component of the conserved neutron (proton) number density current 
(p^). The proton Fermi surface is displaced by K^^. We are interested in the case << kn, kp, but the expressions 
for general are not more complicated than the power series expansion. 
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Noting that 

{i^[j^iid^-g^uj^,)-m,]ij) = (20) 

we find 



,2 ^2 



+ = fc + 5a.u;"^ + ml , (21) 



wliere we have dropped expectation value brackets for the mean values of the fields. The energy e of a baryon in a 
plane-wave state is given by 



s{k) = E{k) - 5„a;° = \/{k + gu,uj'i) +ml- g^uj° . (22) 

Thus we see that wo contributes a constant shift, u)z gives a preferred frame for the momenta, and a renormalizes the 

mass to the Dirac mass. 

As an example of how the expectation values are evaluated, we give the scalar density (letting = K, for ease of 
notation): 

<^^> = jirLf''^ 

' ^ d^k , +A/ d'k , (24) 



(27r)' J|fc|<fc„ ^(fc + 5^a;^?)2+m2 {'^T^f J\k-Kz\<k, ^ + g^^^^)^ + ^2 

Sk , + I d^k , (25) 



and the average four- velocity components of the baryons: 

= I d^k + [ d^k 

(27r) J\k\<k„ (27r) J\k-Kz\<kp 



^ [ d^k + [ d^k , (26) 

Tt) J\k\<k„ (27r) J\k\<k„ 



{2tt) J\k\<k„ {2tt) J\k\<k, 

{2nrJ,,, ' dkz 



k'' + gojOJ' 



-^f d^k , + +A/ ^ 

(27r) i|fe|<fc, ^0 + g^^z^)2^^2 (27r) i|fc-if2|<fc, ^(fc + ^^^^ J)2 + ^2 

2 f 3 fc^+fifo-w^ , 2 f ^3,_ fc^ 5^0;^ + ii- 



3 , d^k + / (i^fc ' ' . (27) 

(27r) J|fe|<fc, ^{h + g^u:^zf+m.f (2^) "^l^l<'=^ {k + g^uj^z + Kzf + ml 

Thus we have reduced the problem to a set of nonlinear equations for the m*, loq, and fields that must be solved 

mimcrically. This can be done for any set of the input parameters A:„, kp, and K. The interaction and mass 
parameters for the effective fields have been determined from nuclear physics, and they are discussed further below. 
Once this is done, we still need expressions for the stress-energy tensor, which is the input for the Einstein equations. 

Again, specializing to the zero-momentum frame of the neutrons, the only non-zero stress-energy tensor components 
are 

m = -\cl {{^I'^f - {i>l'^f) - \c-^ K - ™*) - <^7'fc,V) , (28) 
(TO) = (V^T^fcz^) , (29) 

{TD = (ry) = \cl ((V^lV)' - ^i>l'^f) - \c-^ (m - m.f + <^7"fc.^> , (30) 
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(31) 



Some of the expressions have been simpUfied using the equations of motion. 

Each component of {Tj^) again involves an integration over the Fermi surfaces, but now in terms of completely 
known parameters. For example, to determine (T^^), we need 



(27r)3 
2 

2 

(2^ 



\k\<k„ 



dE 
dk^ 



(27r)3 



d-^kt 



k—Kz\ <kp 



dE 
dk^ 



\k\<k„ 



d^k k^ [kz + guji^z] 



k + g^Lo^z] + ni 



-1/2 



I I <C kp 



d^k [kz + K] [kz + gu^toz + K] 



and for {T'^) 



(27r)3 



k\<k„ 



(27r)3 



\k\<kp 



(k + g^Lu'^z + Kz^ +ml 



fc3 

d^k (kz +K) = -^K 



^1/2 



(32) 



(33) 



The main result of this section is thus a well-defined prescription for producing the functions (T^f) (fc„, fcp, Kz). In 
the next section we take this prescription and produce from it the so-called master function, including entrainment, 
that is used in the general relativistic superfluid field equations. 



III. GENERAL RELATIVISTIC SUPERFLUID FORMALISM 

The formalism to be used here, and motivation for it, has been described in great detail elsewhere [|], H, ^ |l^, 
|l8| , |20| , pT| , and so we will review only the highlights. The central quantity of the superfluid formalism is the 
master function A. It depends on the three scalars = —n^n^, = —p^p^ and = —PfiU^ that can be formed 
from the conserved neutron (ri^) and proton {jp^) number density currents. Furthermore, the master function is such 
that —A{n^,p'^, x"^) corresponds to the total thermodynamic energy density if the neutrons and protons flow together 
(as measured in the comoving frame). Once the master function is provided the stress-energy tensor is given by 

Ti;^^Sii + n^^t,,+p^X. , (34) 

where 

Vp-P^Xp (35) 



is the generalized pressure, and 



^i, = Bn„ + Apu , (36) 
= An„ + Cpu , (37) 



are the chemical potential covectors. We also have 

dk dA OA 

The momentum covectors Hi, and Xu ^.re dynamically, and thermodynamically, conjugate to "ny and p'^ and their 
magnitudes are the chemical potentials of the neutrons and the protons, respectively. The two covectors also make 
manifest the entrainment effect; that is, we see that the momentum of one constituent (/i^, say) carries along some 
of the mass current of the other constituent {^.y is a linear combination of nn and pu)- We can also see that there is 
no entrainment unless the master function depends on x^. 

The field equations for this system take the form of two conservation equations for the neutrons and protons, i.e. 

V,,n^ = , V^p'^^O, (39) 

which is a reasonable approximation given that the weak interaction timescale is much longer than the dynamical 
timescale of neutron stars for small amplitude deviations from equilibrium [p3| , and two Euler equations, i.e. 

^i''V[^A'H=0 ' P^V[,,x.]=0, (40) 



where the square braces means antisymmetrization of the enclosed indices. 
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A. Extracting the Master Function from the Mean Field Results 



The two scales that enter this problem are the microscopic, on the scale of the nucleons, and the mesoscopic where 
one speaks in terms of the two interpenetrating supcrfluids. The fundamental "particles" at the flnid level arc the 
fluid elements which contain, say, an Avogadro's number worth of nucleons. The connection between the micro- and 
meso-scopic levels is via the averaged stress-energy components calculated earlier. Consider a fluid element deep in 
the core of the neutron star and orient the local coordinate frame in such a way that the z-axis of the frame is in the 
same direction as the proton momentum with respect to the neutrons. As shown just below, a unique combination 
of the averaged stress-energy components determined via the mean field theory will yield the master function. As 
this quantity is a scalar, the functional relationship we obtain between A and the two particle number densities and 
the relative velocity of the protons with respect to the neutrons can then be applied anywhere in the star. 

The key idea is to use the (local) relationship 



(41) 



to obtain A. In the perfect fluid case, the identification is made immediate by the fact that there is only one 

four-velocity for the system, and hence a preferred rest-frame for the particles. The local energy density of the 
fluid is thus uniquely obtained from A = — {Tjf) u'^u^. In the superfluid case there are two four- velocities, and 
thus no preferred rest-frame. Fortunately, we can still obtain A in a unique, and covariant, way, by using the trace 
(T) = {T^) and the three scalars that can be formed from contracting (T^) with and p^, i.e. (T^) n^vy , (T^) n^p", 
and {Tlf)pnp'^. We thus find that A is given by 



A = -l{T) + l{x'-nVr' (nV 



p-' 



-x^[n^'p-'+p''n'']){T^.) , 



(42) 



and the generalized pressure is 



In like manner we find that 



* = -((r)-A) . 



A = -{n^p''{Tlf)+x^A)/{x^-nV) , 

B = {p^p''{Ti;)+p^A)/{x^-nV) , 



(43) 



(nX(T^.^)+«'A)/(x4-nV) • 



(44) 



One other necessary component of uniting the mean field theory with the superfiuid formalism is to relate (locally) 
and p^^ to the mean particle fiux of the neutrons and protons; i.e. 



pu<; 



(45) 



where ipn and ijjp are the neutron and proton, respectively, components of the Dirac spinor ip. Recall again that we 
have arranged that the average neutron and proton particle fluxes are in the z-direction. Thus, the unit vectors have 
only two components: 



It thus follows that 

A = (T°) + {T!) - (t;) . 
We will use cylindrical coordinates and define 4>z = Bu^z so that 



(46) 

(47) 



m* = m - -^m* 




dk, [kl + (j)l+ml + 2<j>,k, 



,1/2 



dkz 



k'p + {cP,+Ky+mi + 2{(t>, + K) k. 



1/2 



7 





J 7 

akz 






1 




3^ 




1 


r 


2^ 




1 




2^ 


J ~kp 



r," - —k^ 



1/2 



(k, + (t>z + Ky + ml 



1/2' 



(48) 
(49) 



/ dfc, (fc, + 4>,) [ki + mi+4>i+ 2c^,k 



,1/2 



1/2 



1/2 



/ dk, (fc, + 0, + if ) fc^f^ + mi + (0, + i-T)" + 2 (0, + K) k. 



2 , 2^/2 



(fc, + 0, + if)" + m 



(50) 



It is to be remarked that even though only the protons are given spatial momentum K the neutron four-velocity 
nevertheless has a non-zero spatial component. This, in fact, is a signature of the entrainmcnt effect, which is a 
momentum induced in one of the fluids will cause part of the other fluid to flow. 

Of primary importance to the fluid equations are the A, B, and C coefficients. We could, in principle, use Eq. ( |45| ) 
to express (n^,p^,a;^) in terms of {kn,kp, K), but practically speaking this is not possible. Fortunately we see from 
Eq. (^) that we can construct these coefficients algebraically from the mean field values of the stress-energy tensor 
components. Thus, when it comes to the numerical work, we use the set (fc„,/cp,if) as the independent variables. 
Note that because the master function is a scalar, it must be invariant if if —>■ —K, and is thus an even function of 
K. 



B. Equilibrium models 



The equilibrium configurations are spherically symmetric and static, so the metric can be written in the 
Schwarzschild form 



ds^^-e'^Wdt^ + e^^'^^dr^-^r^ {d9 



2 , -2 

' sm 



The two metric coefficients are determined from two Einstein equations, which are written as 



A' = - — — - Sirre^ A| 



f - e 



X 











(51) 



(52) 



r r 
The equations that determine the radial profiles of n{r) and p{r) have been derived by Comer et al. pi| and they are 



"^o|oP'+^S|o"' + ^mIo^' = , CS|,p'-|-^o|„n' + ixlo^' = 



where 



Bl ^ B- 



op'' 

^dB 2 
on'' 

2— 2 4 



din? 



2^n" + 2— p 



dA 2 
dp^^ 



dA 



pn 



On'' 



dA 2 
dA , 



^dA^ ^ 

9p2 ^ dx^ 



(53) 

(54) 
(55) 
(56) 



and A|q, ^|g, /ilp, and x\o &re given in the appendix. A zero subscript means that after the partial derivatives are 
taken, then one takes the limit if ^ 0. 

Of course, since the variables that are more suited to the mean field theory are the two Fermi wavenumbers, kn 
and fcp, we replace everywhere n = k'^/iir^ and p = fc^/STr^ and solve for the wavenumbers instead. We have also 
found a more convenient way of determining the Dirac effective mass |q (fc„,/cp) and that is we have turned the 
transcendental algebraic relation in Eq. (pj) of the appendix into a differential equation via 



m 



*io 



dm^, 
dkn 



k^ 



dm^, 
dk„ 



kp 



(57) 
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where k'^ and k'^ are obtained from Eq. (|5^). 

The "boundary" conditions that must be imposed include a set at the center and another at the surface of the star. 
Demanding a non-singular behavior at the center of the star imposes that A(0) = 0, and consequently that A'(0) and 
i^'iO) must also vanish. This and Eq. (|5^ ) imply that fc^(O) and fcp(O) have to vanish as well. A smooth joining of 
the interior spacetime to a Schwarzschild vacuum exterior at the surface of the star, i.e. r = R, implies that the total 
mass M of the system is given by 

M = -All i Ar [r) (58) 

Jo 

and that (i?) = 0. 



C. The low velocity limit for fluid elements 

Two immediate applications of the formalism developed here are to model slowly rotating configurations 0] and 
linearized perturbations, or quasinormal modes |^. In both cases the fluid element velocities are small in the 
sense that they are typically only a few percent of the speed of light. The net effect is that these applications require 
only the first few terms from an expansion of the master function in terms of the entrainment parameter (i.e. in 
the canonical formulation, and in what we present here). Such an expansion has been described in |Q, ||, and 
thus only the highlights will be reproduced here. It should be noted, however, that if one wanted to model rapidly 
rotating superfiuid neutron stars |24|, say, then the expansion to be described below will be inappropriate. 

For a region within the fluid that is small enough that the gravitational field does not change appreciably across it 
one can show that 

X ^ np\ , , . (59) 

If it is the case that the individual three-velocities Vn,p are small with respect to the speed of light, i.e. 

^ « 1 , (60) 

c 

then it will be true that ~ np to leading order in the ratios it„/c and Vp/c. Thus, an appropriate expansion of the 
equation of state is 

oo 

A(n2,p2^a:2)^^A,(n2,p2)(^2_^p)^ ^ (g^) 

i=0 

since x^ — np is small with respect to np. In this case the A, ^q, etc. coefficients that appear in the field equations 
can be written as 

oo 

A = -'^i \i{n^ ,p^) (yX^ - npY , 
1=1 

n on n n ^ — ^ on 

i—l 

C = -'-^-^A-'-±^{x^-npy , 
p op p p op 



opon opon ^ ' 

i—l 
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For quasinormal mode and slow-rotation calculations, each of the coefficients are evaluated on the background, so 
that — np, and thus only the first few are needed. In fact, one needs to retain only Aq and Ai, where the latter 
contains the information concerning the entrainment effect. Some details are given in the appendix, and the final 
results are 



CL = 



«0| 



37r2fc5 k 

XIq _ 2 ^ _ 'ii^n 
1.3 1,3 K ,,2 1 1,3 




, (63) 
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5 /x^Iq kp 



kp 




37r2 



1,2 1,3 



iL2p 



1,2 
ftp 



TT 



2,5 



dk 
dk,i 

dA 
dk„ 



d^A 

" dki 



^2 kn + m*|g 



lo 9fc„ 



)L2 
ftp 



fc2 -t- TO2|g 



(64) 



(65) 

(66) 
(67) 
(68) 



where dm^/dkn\Q and dm^/dkp\^ can be found in the appendix and rrig = m/1836. For reasons to be discussed 
below, we have included contributions due to a normal fluid of highly degenerate electrons. 



D. Equilibrium configurations 

We now use our model to construct static and spherically symmetric configurations. A priori there are two input 
parameters, which are the neutron and proton Fermi wavenumbers at the center of the star. However, we can reduce 
this to just the neutron wavenumber by imposing the condition of chemical equilibrium between the nucleons. In 
order to have a chemical equilibrium that is believed to be representative of neutron stars (i.e. proton fractions 
= p/i'iT' + p) ~ 0.1), we have added to the master function a term (see, for instance, ||, ^) that accounts for a 
highly degenerate gas of relativistic leptons (in our case, just electrons). Fig. gives the mass M as a function of 
the central neutron number density n(0). We see the typical behaviour of general relativistic neutron stars, and that 
is a maximum value for the mass. Beyond this maximum, the stars will be in unstable equilibria. As canonical 
models of superfluid neutron stars, we have chosen configurations that are near to the maximum mass, but on the 
stable branch of the curves (cf. Table |). 

In several earlier studies [Q, |^, a parameterized model for entrainment has been used that is based on the 
Newtonian calculations of Borumand et al. [Tl|| . The parameter, £mom, in this model can be shown to be related to 
our coefficient via 

-^'"^ (69) 



m — {n +p) 



In the previous studies the "physical" range has been taken to be 0.04 < emom < 0.2. But already with this formula 
we find that emom has a singularity whenever m = A\q (n + p). And in Fig. || we see that indeed Smom has a 
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FIG. 1: Mass M (in units of solar mass Mq) vs the central neutron number density n(0) (in units of fm for the mean field 
coupling values of Models I and II of Table 



TABLE I: Parameters describing our choice of mean field and canonical superfluid neutron star models. The two values for 
and represent the two extremes given in that have been determined from nuclear physics. Note that the baryon mass 
is m = 4.7582 fm"\ 



Model I 


Model II 


cl = 12.684 
cl = 7.148 
1^(0) = -2.316408 
fc„(0) = 2.8 fm"^ 
Xp{0) = 0.101 
M = 2.509 Mq 
R = 11.696 km 


4 = 8.403 
cl = 4.233 
1^(0) = -2.288385 
fc„(0) = 3.25 fm-i 
Xp{0) = 0.102 
M = 1.996 Mq 
R = 9.432 km 



discontinuity in both curves (which are for Models I and II of Table |). This is perhaps the most significant short 
term result of this work, and hence taking £mom to be a constant must be considered a non-viable option. 

An alternative parameterization is that of Prix et al. . They point out that there is some ambiguity in what is 
meant by the nucleon effective masses (i.e. the Landau, as opposed to the Dirac, effective masses |l^), which can be 
traced to whether one chooses to define these masses with respect to the zero-momentum or zero-velocity frame of 
the nucleons. Their parameterization makes use of the zero-velocity frame. In this case one finds 



A\q n 
m 



(70) 



Immediately we see that this formulation does not have an obvious singularity anywhere, and this is reflected in Fig. ^ 
which contains plots of the radial profile of Syei for Models I and II of Table |. The "physical" range has been taken 
to be 0.4 < Eyf^i < 0.7. We see from the figure that this is the range in the outer portion of the star, but in the core 
the entrainment stretches well out of this range. 

Also used in Prix et al. H is the so-called symmetry energy parameter 



d^i/dp 
dx/dn 



(71) 



that can be related to terms in the equation of state that tend to force an equal number of protons and neutrons 
(as in most nuclei). For the relativistic mean field model used here, it is not difficult to show that a = \ (which is 
consistent with the range of values used by 0). 
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FIG. 2: The entrainment parameter Smom as a function of radius for Models I (left panel) and II (right panel) of Table ^. 
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FIG. 3: The entrainment parameter £„ei as a function of radius for Models I (left panel) and II (right panel) of Table | 
IV. SLOW ROTATION CONFIGURATIONS: THE FRAME-DRAGGING 



The key distinguishing feature of slowly rotating superfluid neutron stars is that the neutrons can rotate at rates 
different from that of the protons. The slow-rotation approximation is valid when the angular velocities are small 
enough that the fractional changes in pressure, energy density, and gravitational field due to the rotation are all 
relatively small. This translates into the inequalities (cf. ||l|, |27| ) 



nl or or n^np « ^ , 



(72) 



where the speed of light c and Newton's constant G have been restored, and R and M are the radius and mass, 
respectively, of the non-rotating configuration. Since GM/c^R < 1, we also see that 



ftnR « c and flpR << c , 



(73) 



and thus the slow-rotation approximation ought to be useful for most astrophysical neutron stars. In fact a comparison 
[Q, 1^] of the above conditions to empirical estimates for the Kepler frequency (i.e. the rotation rate at which mass- 



12 



shedding sets in at the equator) that can been obtained from calculations using realistic supranuclear equations of 
state reveals that even the fastest observed pulsars can be classified as slowly rotating. 

The only quantities that contain terms linear in the angular velocities are the metric coefficient uj{r), that represents 
the dragging of inertial frames, and the fluid four- velocities. All other effects due to rotation enter at the second-order 
in the angular velocities. It is useful to define 



LU - Cl„ 



Lp — CO — rip 



(74) 



Up to an overall minus sign, these represent rotation frequencies as perceived by local zero-angular momentum 
observers. The Einstein equation that determines the frame-dragging has been shown to be ||l| 



(75) 



It is of the same form as that obtained by Hartle |^ except for the non-zero source term on the right-hand-side. 

Exterior to the star, there is vacuum, and so the solution for the frame-dragging is the same as that considered by 
Hartle |27l, i.e. 



2J 

Assuming that the frame-dragging is continuous at the surface of the star, then 



(76) 



J 



Stt 



(77) 



where J is the total angular momentum. Andersson and Comer Q have furthermore shown that the neutron total 
angular momentum is 



J n 



Stt 



d^^4g(A-.)/2 



+ A\ 



(78) 



and 



Jp — 



Stt 



xlaPLp + Al^np (rip - fin) 



(79) 



for the proton total angular momentum, from which it follows that 



(80) 



Any solution of Eq. (|75|) for the frame-dragging is to be such that the interior matches smoothly onto the known 
vacuum solution in Eq. (|7q). This means that we must have, for instance, 



Lp{R) 



2J 



(81) 



We can easily see that Lp and its derivative are smooth provided that we have 



Lp{R) 



R dLp 



(82) 



r=R 



Having obtained a value for Lp{0) that satisfies Eq. ( |8^), a n acceptable solution is in hand, and we can thus determine 
the angular momentum of the configuration from Eq. ( pi|) . 

In Fig. 1^ we have plots of the radial profile of the frame-dragging for Model I for a range of values of the ratio 
n„/J7p. For the values considered in the left panel we see that the frame-dragging is much like that of an ordinary 
one-fluid star, and is consistent with solutions obtained by Andersson and Comer For the negative ratios, we 
see that the frame-dragging is negative but increases monotonically towards zero. This is the behaviour we should 
expect, since the bulk of the matter is simply rotating the opposite way. There is some asymmetry between the 
negative and positive ratios, but that is due to the small number of protons that rotate oppositely to the neutrons 
when the ratio is negative. In the right panel, we examine the solutions near to a ratio of zero. The frame-dragging 
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FIG. 4: The radial profile of the frame-dragging Ld{r) for Model I of Table |. In the left-panel we have curves for Q,„/Q,p = 
(-2.0, -1.0, 0.0, 1.0, 2.0), and in the right we have taken Q^/Qp = (-0.11, -0.10, -0.09). 



is no longer monotonic and actually becomes negative inside the star. An explanation of this can be understood as 
follows: in the interior the protons carry most of the angular momentum and thus have the largest impact on the 
frame-dragging, but further away from the center, the much larger mass contained in the neutrons begins to dominate 




Fig. |5| considers the same range for the ratio of the angular speeds, by showing how the total angular momentum 
J, and the neutron and proton angular momenta, Jn and Jp respectively, vary as Qn/^p is changed. As one might 
expect, when the ratio becomes greater than one, then the angular momentum in the neutrons is significantly greater 
than that of the protons. Likewise, as the ratio becomes smaller we find that the protons can dominate. Something 
of a surprise is the extreme right-hand-side of the curves, where Jp actually becomes negative, and yet the angular 
speeds of the neutrons and protons have the same sign. We can explain this behaviour as a purely general relativistic 
effect that is intimately connected with the frame-dragging. With respect to infinity, particles that are rotating 
around at the same rate as the local inertial frames are found to have zero angular momentum. Thus, those particles 
that would be lagging behind the frames, even though their angular trajectories would be in the same direction as 
the frames, will nevertheless have negative angular momentum. Finally, one other feature is the configuration where 
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FIG. 5: The neutron J„, proton Jp, and total J angular momenta vs the ratio fin/fip for Model I of Table |. 

Jn = Jp- In this case, the angular speeds are not equal, nor are the total neutron and proton particle numbers equal, 
and yet the angular momenta of both fluids are the same. 

V. CONCLUSIONS 

We have developed a formalism that uses relativistic mean field theory for supra-nuclear density matter that can 
be applied to general relativistic superfluid neutron stars. In this formalism we have also allowed for the entrainment 
effect between the various superfluid species. We have also shown how to use our formalism in the relativistic 
superfluid field equations that have been recently developed for modelling slowly rotating equilbirium configurations 
[Q, as well as the linearized oscillations Our results indicate that parameterized models of entrainment that 

use a zero-momentum rest-frame definition for nucleon effective masses have a singularity, whereas those based on 
zero-velocity rest-frames do not. 

Our results should find a wide range of applications, not the least of which is to understand better the role of 
entrainment in the superfluid modes of oscillation (e.g. the avoided crossings described by [^) and subsequent imprints 

^ that may be left in neutron star gravitational waves (emitted, for instance, during glitches). Applications planned 
for the near future will include numerical studies of rapidly rotating superfluid neutron stars (using an adaptation of 
the very accurate LORENE code [^) and continued research on the newly discovered two-stream instability ||^, |j, 
which has been proposed as a trigger mechanism for glitches in pulsars. For the rapid rotation calculations one 
must necessarily imploy the full formalism discussed here in the sense that K will no longer be kept small, since 
the LORENE code is specifically designed to accurately handle relative velocities of the neutrons with respect to the 
protons that approach the speed of light. And for the two-stream instability entrainment provides one of the main 
couplings between the two fluids. 
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Appendix: Limiting Forms 



The slow-rotation approximation is such that only terms up to and including 0{fl^ p) are required. This translates 
into keeping only those terms in the mean field theory up to and including 0{K^). This is because those quantities 
like m* and A are scalars, and can only depend on terms that are even in K. Likewise, those quantities that are 
like vectors, e.g. Wz, can only depend on terms that are odd in K. Because m* and lOz are known only implicitly, we 
determine their expansion coefhcients by assuming they take the form 



dK 



K 



where 



m — m* 



27r2 



1.2 



+ rnl\g + kpJk^+ m2|p + - m^l^ln 



— kn + \/kn + ™*lo 

kn + \/k\ + Iq 



lo 



In 



'*io 



fc2 + m% 



By inserting Eq. ( p3[ ) into Eq. (pO]), and expanding and keeping terms to the appropriate orders, we find 

r2 fc3 kl \^ 

dK - " ^ 



dkn 
9m, 



dkn 




3^^' ^/kl+ml\, dK 
1 



l3 



37r2 



fc2 + m 



2 

* lo 



dK 



1 K 



(83) 



(84) 



(85) 



(86) 



(87) 



(89) 



We note that the coefficient d'm,,/dK^\^ cancels everywhere, which is why it is not written here. Also, we find 
A|o = 




Xlo = 9^^f) + \/kp + ml\^ , 



(90) 

(91) 
(92) 
(93) 



16 



The condition of chemical equihbrium is that /^lo = xlo + \/^'p 
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